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Abstract 

We address the problem of the relative importance of the intrinsic chaos and the 
external noise in determining the complexity of population dynamics. We use a re- 
cently proposed method for studying the complexity of nonlinear random dynamical 
systems. The new measure of complexity is defined in terms of the average number 
of bits per time-unit necessary to specify the sequence generated by the system. 
This measure coincides with the rate of divergence of nearby trajectories under two 
different realizations of the noise. In particular, we show that the complexity of a 
nonlinear time-series model constructed from sheep populations comes completely 
from the environmental variations. However, in other situations, intrinsic chaos can 
be the crucial factor. This method can be applied to many other systems in biology 
and physics. 
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1 Introduction 



Recently several outstanding papers[l,2,3,4,5,6,7,8] have applied physical and 
mathematical methods to ecology and population dynamics. This is a very 
important development. In fact, interdisciplinary research can produce very 
significant ideas. 

There exists a great controversy in ecology[9,10,ll,12,13,14,15] concerning the 
relative importance of intrinsic factors and external environmental variations 
in determining populations fluctuations. In this article we address this prob- 
lem using a recently proposed method[16,17] for studying the complexity of a 
nonlinear random dynamical system. This method characterizes the complex- 
ity by considering the rate K of divergence of nearby orbits evolving under 
two different noise realizations. We can show that this measure is very effective 
for investigating nonlinear random systems. In Ref . [14] a nonlinear time-series 
model is constructed from sheep populations on two islands in the St. Kilda 
archipelago [18, 19]. We investigate the complexity of this model using the new 
technique. We have shown that the complexity of the system comes com- 
pletely from the environmental variations. This combination of new methods 
is a very powerful tool for quantifying the impact of environmental variations 
on population dynamics and can be applied to other systems. 

The paper is organized as follows: In Section 2 we recall the definition of 
complexity for random dynamical systems. In Section 3 we recall the definition 
and properties of the random sequences given in Ref. [20,21,22,23], and we 
compute the complexity for the random sequences and a particular random 
map. In Section 4 we discuss a nonlinear time-series model constructed from 
sheep population data and we compute the complexity for this model. In this 
section we also show that for a generalized model, the complexity can depend 
on both, the intrinsic chaos and the environmental variations. In Section 5 we 
briefly discuss some aspects of the problem of distinction between deterministic 
chaos and noise. 



2 Complexity in random dynamical systems 

Recently a new measure of complexity was introduced[16,17] in terms of the 
average number of bits per time-unit necessary to specify the sequence gen- 
erated by the system. This definition becomes crucial in random nonlinear 
dynamical systems as the following 



X n+ i — f(X n , I n ), 
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where I n is a random variable (e.g. noise). This measure coincides with the 
rate K of divergence of nearby trajectories under two different realizations of 
the noise. The method of calculating the Kolmogorov-Sinai entropy with the 
separation of two nearby trajectories with the same realization of the noise 
can lead to incorrect results[16,17]. The complexity of the dynamics generated 
by (1) can be calculated as 

K = A0(A) + h, (2) 

where A is the Lyapunov exponent of the map, which is defined as 

A= lim vT 1 In \Z n \, (3) 



where Z n+1 = (df(X n )/dX n )Z n , h is the complexity of I n , (which can be cal- 
culated as the Shannon entropy of the sequence I n ), and 6(X) is the Heaviside 
step function. This function is defined as follows: 6(X) = if A < 0; 6(X) = 1 
if A > 0. For a detailed explanation of the relationship between the definition 
of complexity as the average number of bits per time unit necessary to specify 
the sequence, the rate of divergence of nearby trajectories under two differ- 
ent realizations of the noise and the equation (2), see Ref.[36]. The definition 
in this form was given in the original papers [16, 17]. However, in a different 
formalism Eq.(2) could be considered as the starting definition of K. On the 
other hand, there are many alternative measures of complexity. So we should 
check the effectiveness of this new method. In the next section, using some 
random sequences and a random map, we will show that the rate of divergence 
of nearby trajectories under two different realizations of the noise indeed can 
be calculated using equation (2). 



3 Random sequences 



Very recently[20, 21,22,23] we have investigated explicit functions which can 
produce truly random numbers 

X n = sin 2 \6tx Z n ). (4) 



When Z is an integer, function (4) is the exact solution to chaotic maps. How- 
ever, when Z is a generic fractionary number, this is a random function whose 
values are completely independent. Using these functions (or an orthogonal 
set of them) we can find exact solutions to random maps as equation (1). 

Let us discuss some properties of function (4). Let Z be a rational number 
expressed as Z = p/q, where p and q are relative prime numbers. We are 
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going to show that if we have m + 1 numbers generated by function (4): 
X , Xi, X2, X 3 , X m (m can be as large as we wish), then the next value 
X m+ i, is still unpredictable. This is valid for any string of m + 1 numbers. Let 
us define the following family of sequences 



X n k ^ = sin 2 [tt (6 + kq m ) (p/q) n ] , 



(5) 



where k and m are integer. For all sequences parametrized by k, the first m + 1 



values are the same. This is so because Xi k ' m ^ = sin 2 



n6 Q {j)/q) n + nkp n q 



n„(m—n) 



sin 2 [7i9o(p/q) n ], for all n < to. Nevertheless, the next value 



X 



(k,m) 
m+1 



= sin 



K0 o (p/q) 



m+1 



+ 



nkp 



m+1 



(6) 



is unpredictable. In general, X^+1 can take q different values. These q values 
can be as different as 0, 1/5, 1/2, V2/2 or I. For Z irrational there can be an 
infinite number of different outcomes. 

Function (4) with Z = 3/2 (i.e. X n = sin 2 [#7r(3/2) n ]) is a solution of the 
following map 



X n +i — — 



l + / n (l-4X„)(l-X n ) 



1/2 



(7) 



where 



In = 



cos [9n (3/2) r 



{l-sin 2 [^7r(3/2) n ]} 



1/2 



(8) 



if sin 2 [9n (3/2)"] ^ 1; and 



/„ = !, 



(9) 



if sin 2 [dix (3/2)"] = 1. 

A careful analysis of function I n yields that I n is an unpredictable function 
that takes the values ±1 with equal probability. A particular realization of 
I n for 9 = 0.77 is the following: 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, -1, -1, 1, 
1, -1, 1, -1, 1, ... The same analysis of Eq.(5) made above (in this case for 
Z = 3/2) confirms these results. On the other hand, a statistical investigation 
of the outcomes of function I n corroborates these findings. In fact, it does not 
matter how many past values I , I±, I m we already know, the next value 
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cannot be determined. It can be either 1 or -1 with the same probability. In 
other words, I n behaves as a random coin toss. 

Now we can check some of the results discussed by the authors of Ref.[16,17]. 
In the case of the random map (7) A = ln(3/2) and h = In 2. Thus, K = In 3. 
Here we give a brief explanation of these results. The map (7) can be re-written 
on the following form 



X n +i — 



l + (l-4X„)(l-X„) 1 /2 
1-(1-4X„)(1-X„) 1 / 2 



with probability |, 
with probability \. 



(10) 



This is a form also compatible with the application of equation (2) (See 
[16,17]). After the transformation X n = sin 2 (F n ), both equations 



X 



n+l 



l + (l-4X n )(l-X n ) 1/2 



'II) 



and 



n+l 



l-(l-4X n )(l-X n ) 1 / 2 



(12) 



can be converted into piece-wise linear maps where the absolute value of the 
slope \dY n+1 /dY n \ is constant and equal to 3/2. 

Using the equation (3) for the Lyapunov exponent, we obtain the exact value 
A = ln(3/2). The Lyapunov exponent is invariant with respect to the transfor- 
mation X n = sin 2 (y„). If we calculate numerically the Lyapunov exponent of 
maps (7), (11) and (12), we also obtain the value A = ln(3/2) approximately. 

Considering the properties of the sequence I n (that takes the values 1 and —1 
with equal probability), it is trivial to get that h = In 2. Applying equation 
(2), we obtain K = In 3. All these calculations have been made using the 
definitions of the quantities and the algebraic structure of equation (7). Now 
let us consider the analytical solution of map (7): 

X n = sin 2 [07r(3/2) n ] . (13) 



If we investigate equation (13), it is possible to prove that, on average, for a 
given 9, nearby trajectories will separate following the law d ~ (3/2) n , where 
d is the distance between the trajectories. This yields that A = ln(3/2), which 
corroborates a previous result. 



5 



A more important calculation is that of K. We wish to compute the rate of 
divergence of nearby trajectories under two different realization of the noise. 
In the "language" of the exact solution (13) this is equivalent to investigate 
the average divergence of trajectories that are very close for n = 0, but with 
different values of 9 (recall that different realizations of the random variable 
/„ (equation (8)) are produced with different values of 9. This analysis yields 
the following result K = In 3. And this is a corroboration of equation (2) for 
this system. 

Now let us resort to numerical calculations. We have produced numerically 
10000 values of X n using both, the dynamical system (7) and the function (13). 
Then, we have computed the complexity of these sequences using the Wolf's 
algorithm [24]. The result is very close to In 3 (in fact K 1.098). Moreover, 
even an independent calculation of the complexity of this dynamics using 
different methods[24,25] produces the same result. In Ref. [25] a new method 
for the calculation of the complexity of a sequence is developed. This method 
has been shown to be very effective for the calculation of complexity of finite 
sequences[25, 26, 27,28]. We start with a sequence of values Lq, U2, U3, Un] 
from which we can form a sequence of vectors 

X(i) = [U i ,U i+1 ,...,U i+m - 1 ]. (14) 



Now, we will define some variables: 

nmt \ _ ( num ber of j such that d [X(i),X(j)] < r) 

6 * (r) " (N-m+1) ' (15) 



where d[X (i) , X (j)] is the distance between two vectors, which is defined as 
follows: 

d[X(i),X(j)] = max(\U i+k -! - (k = 1,2, ..,m). (16) 



Another important quantity is 

In C™(r) 
N-m+1 



N-m+l 1 rim/ r \ 

r(r)= e - l ( ;v (17) 

1=1 



Using all these definitions, we can calculate the complexity 

K{m,r,N) = <p m {r) -0 m+1 (r). (18) 



This measure depends on the resolution parameter r and the embedding param- 
eter m, and represents a computable framework for the "Shannon's entropy" 
of a finite real sequence. 
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It is interesting that when we calculate numerically the complexity of the 
sequences produced by the map (7) and the exact solution (13) (with r = 
0.025, and different m > 2) we obtain K ps 1.098. 

Using functions (4) (with different values of Z) we can also solve maps where 
the Lyapunov exponent is negative and, nevertheless, due to the existence of 
external noise, the complexity is positive. In the presence of random pertur- 
bations, K can be very different from the standard Lyapunov exponent and, 
hence, from the Kolmogorov-Sinai entropy computed with the same realiza- 
tion of the noise. 

In general, if we apply the measure of complexity K to our functions (4), then 
we obtain the following results: for Z = p/q, K = In p. If Z is irrational, the 
complexity is infinite. We should add some comments about these computa- 
tions. When Z is integer, function (4) is equivalent to a univalued chaotic 
map of type X n+1 = f(X n ). In this case, this measure coincides with the 
Kolmogorov-Sinai entropy i.e. K = A. So A = InZ. When Z = p/q, where 
p and q are relative primes, function (4) produces multivalued firsts-return 
maps[20,21,22,23]. In this case, the information lacking is not given by K = A, 
but is larger: one loses information not only in each iteration due to A > 0. 
One has also to specify the branch of the map (X n , X n+ i) in each iteration. 
We should apply the formula (2), where h is the entropy of the random jumps 
between the different branches of the map (X n ,X n+ i). For Z = p/q, there 
are q branches in the map (X n , X n+ i). Investigating the properties of function 
(4) we arrive at the conclusion that all the branches possess the same prob- 
ability in the process of jumping. This leads to the equality h = Inq. Thus, 
K = \n(p/q) + lng = lnp. The complexity can be obtained by computing the 
separation rate of nearby trajectories evolving in two different realizations of 
the noise. In the case of sequences produced by function (14), this is equivalent 
to using two different 9 for which (at n — 0) the trajectories are close. Such 
procedure exactly corresponds to what happens when experimental data are 
analyzed with the Wolf et al. algorithm [24]. When we apply the Wolf et al. 
algorithm to the sequences generated by our functions and the mentioned dy- 
namical systems, we obtain the expected theoretical results. The complexity 
of the sequences produced by function (4) also can be calculated using random 
dynamical systems for which function (4) is the exact solution as we did in 
the case Z = 3/2. For different values of Z, the result is again K = lnp. 



4 Sheep population model 

In a beautiful work Grenfell et al. [14] used the unusual situation of time series 
from two sheep populations that were very close (and so, they shared approx- 
imately the same environmental variation as for example rain, temperature, 
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etc) but which were isolated from each other, i.e. these populations did not 
interact, in order to study the interaction between noise and nonlinear pop- 
ulation dynamics. They found high correlations in the two sheep populations 
on two islands in the St. Kilda archipelago. They were able to express X n+1 
as a function of the previous population size: X n+ i = f(X n ) + e n +i, where e n 
represents the noise, which is related to the environmental variables (n is the 
discrete time). Here X n = logN n , where N n is the population number. They 
fit a nonlinear self-exciting threshold autoregressive (SETAR) model[29,30,31] 
to the Hirta island (one of the islands of the St. Kilda archipelago) time series. 
The best-fit model is 



X n +\ — &o + b X n + e^+i, X n < c, 



— "1 T 



X n > C, 



(19) 



where c = 7.066; a = 0.848; b = 0.912; a = 0.183; ai = 7.01; a x = 0.293. 
Here a"o,i is the variance of e n . The noise e n is defined as a sequence of in- 
dependent and identically distributed normal random numbers with mean 
and variance a. The model captures the essential features of the time-series, 
including the map X n+ i versus X n . Now we apply the measure of complexity 
K (Eq.(2)) to the model given by equation (19). It is straightforward to show 
that A < 0. This can be done even analytically using equation (3). Thus, the 
complexity of the dynamical system (19) is K — h, where h is the complexity 
of the noise e n . That is, all the complexity of this dynamical system comes 
completely from the environmental variations. We could say that, in this case, 
the extrinsic environmental variations are much more important than the in- 
trinsic factors in determining population size fluctuations. 

This result confirms the results of Ref . [14] . However, we should note that their 
research is based on the very particular situation where we have synchroniza- 
tion of two population fluctuations at separate, but not too distant locations. 

Here we should explain briefly how Grenfell et al. [14] obtained their results. 
They found that the fluctuations in the sizes of the two populations are re- 
markably synchronized over a 40-year period. They explain this synchroniza- 
tion using the fact that the two populations are exposed simultaneously to the 
same environmental variations. Assuming that the same model applies to both 
islands, they use it to estimate the level of correlation in environmental noise 
required to generate the observed synchrony in population fluctuations. They 
found that very high levels of noise correlation are needed to generate the 
observed correlation between the sheep populations on the two islands. They 
also studied observed large-scale meteorological covariates like monthly wind, 
rain and temperature. From this analysis they conclude that the extrinsic 
influences are very important in this particular case of population dynamics. 
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On the other hand, our method can be applied to any other population dy- 
namics. Even if we have only one isolated population in the same region. 
The research program is the following: the data should be fitted by a SETAR 
model and, after that, the complexity can be calculated using equation (2). 
In fact, many nonlinear population models can be approximated by a SETAR 
model. This is a very clear theoretical result. It does not depend on further 
statistical assumptions or approximate investigations. Once we have recon- 
structed the model from the data (and this is a step that we cannot avoid 
in any other method), we can prove rigorously that the Lyapunov exponent 
is negative and that all the complexity comes from the external random per- 
turbations. Our results explain why the environmental variations are more 
important in this particular case. This is due to the density-dependent re- 
lationship X n+1 = f(X n ). In fact, for other animal populations the best-fit 
model can be very different. The form of the density dependence is crucial. 
For instance, suppose that the best-fit model is similar to that presented in 
Ref.[15]: 



X 



n+1 



r + X n + e£k, X n < c, ^ 

( r + 6c ) + (l_ & )X„ + ei 2 | 1 , X n >c. 



Here we should add a short explanation of the origin of Eq.(ll). In Ref. [32] 
Mayrand Smith and Slatkin discuss the so-called MSS model 

N = RNn (21) 



where N n is the population size at time n, R is the maximal net population 
growth rate, b is a measure of the strength of the density dependent reduction 
of the net population growth rate, and N c is the carrying capacity (that is, 
the maximum population size that can be sustained by the area under study). 
If we introduce the transformation X n = lnN n (and the effect of noise), we 
can obtain equation (20) as an approximation, where r = In R, c = In N c (See 
Ref. [15]). However, in the same way as Eq.(19), Eq.(20) can be obtained as a 
SETAR model reconstruction using a given time-series of the evolution of cer- 
tain animal population. In general, as pointed out by Stenseth and Chan[15], 
many nonlinear population models may, on the log-scale, be approximated by 
a dynamical system similar to equation (20). 

For large values of parameter b, the Lyapunov exponent can be positive. In this 
case, both the intrinsic and the external factors contribute to the variability 
of the dynamics. Moreover, it can be that the intrinsic chaotic factors are the 
most important in determining population fluctuations. Nevertheless, we have 
shown that randomness is crucial in ecological models. 
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5 Chaos and noise 



In the research presented in this paper the question of chaos or noise is very 
relevant both in the problems related to the random functions and in the 
problems of characterizing experimental time-series. 

We should say here that recently several important papers have been dedi- 
cated to the question of distinguishing between generic deterministic chaos 
and noise[33,34,35,36,37]. Several limitations have been found for the usual 
methods that are based on the calculation of the Lyapunov exponent and the 
Kolmogorov-Sinai entropy [37]. Many of the practical problems are related to 
the fact that these quantities are defined as infinite time averages taken in the 
limit of arbitrary fine resolution. 

New very strong methods have been developed based on different concepts. 
Some of these methods[33,34,35] are based on the differences in the predictabil- 
ity when time-series is analyzed using prediction algorithms. 

In Ref. [37] this problem is solved by introducing the (e, r) entropy (h(e, r)), 
which is a generalization of the Kolmogorov-Sinai entropy with finite resolu- 
tion e, and where the time is discretized by using a time interval r. 

If the Kolmogorov-Sinai entropy can be calculated exactly and is finite, then 
we can assure that the time-series was generated by a deterministic law. Usu- 
ally the (e, r)-entropy displays different behaviors as the resolution is varied. 
According to these different behaviors one can distinguish deterministic and 
stochastic dynamics. We can even define a certain range of scales for these 
phenomena. 

For a time-series long enough, the entropy can show a saturation range. For 
e — > 0, one observes the following behaviors: h(e) ~ const for a deterministic 
system, whereas h(e) ~ — In e for a stochastic system. In general, predictabil- 
ity can be considered as a fundamental way to characterize complex dynamical 
systems[36]. We have used a method developed in the papers[33,34,35], in or- 
der to investigate numerically the randomness of functions (4). This technique 
is very powerful in distinguishing chaos from random time series. The idea of 
the method is the following. One can make short-term predictions that are 
based on a library of past patterns in a time series. By comparing the pre- 
dicted and actual values, one can make distinctions between random sequences 
and deterministic chaos. For chaotic (but correlated) time series, the accuracy 
of the nonlinear forecast falls off with increasing prediction-time interval. On 
the other hand, for truly random sequences, the forecasting accuracy is inde- 
pendent of the prediction interval. If the sequence values are correlated, their 
future values may approximately be predicted from the behavior of past values 
that are similar to those of the present. For uncorrelated random sequences 
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the error remains constant. The prediction accuracy is measured by the coef- 
ficient of correlation between predicted and observed values. For deterministic 
chaotic sequences this coefficient falls as predictions extend into the future. 
Suppose we have a sequence u±, u 2 , ujy. Now we construct a map with the 
dependence of u n (predicted) as a "function" of u n (observed). If we have a 
deterministic chaotic sequence, this dependence is almost a straight line, i.e. 
u n (predicted) ~ u n (observed) (when the forecasting method is applied for one 
time step into the future). When we increase the number of time steps into the 
future, this relation becomes worse. The decrease with time of the correlation 
coefficient between predicted and actual values has been used to calculate the 
largest Lyapunov exponent of a time series[34]. We have applied this method 
of investigation to our functions(4)[20]. When Z is an integer (Z > 0), the 
method shows that the function (4) behaves as a deterministic chaotic system. 
If Z is irrational, the correlation coefficient is independent of the prediction 
time. Even when the method is applied with prediction time interval m — 1, 
the correlation coefficient is zero (the map (u n (predicted),u n (observed)) cov- 
ers completely the square 0<x<l,0<y<l, showing no patterns). This 
shows that the corresponding time series behaves as a random sequence. When 
Z = p/q, function (4) behaves as a system with both, deterministic chaos and 
noise. In this case, it is better to complement the study with several alternative 
methods. 

As pointed out by Cencini et al. ,[37], all these methods have in common that 
one has to choose certain length scale e and a particular embedding dimen- 
sion m. Thus the scenarios discussed in [36,37] can be very useful in all the 
investigations aimed at the distinction between chaos and noise. 



6 Conclusion 

Cohen [38] has reported that the solutions of chaotic ecological models have 
power spectra with increasing amplitudes at higher frequencies. This is in con- 
trast with the spectra presented in natural populations which are dominated 
by low-frequency fluctuations. Some authors [11] suggest that this is a man- 
ifestation of the interaction between biotic factors and climatic factors. This 
problem shows the difficulties in deciding whether natural populations fluc- 
tuations are determined by internal biological mechanisms or they are mostly 
the result of external environmental forcing. We think that our results can 
help to shed light on this issue. Recently there have been reports [39, 40] on 
population dynamics where the variability of the population originates from 
both deterministic chaos and stochastic processes. The complexity given by 
equation (2) can help to determine the relative weight of both factors. In 
fact, the understanding of the interaction of both deterministic and stochastic 
processes is crucial to model correctly the dynamics of an ecological system. 
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We propose a combined approach to this issue: the SETAR model and the 
new method for calculating complexity. In the particular case of the sheep 
populations in the St. Kilda archipelago, it seems that the population fluctua- 
tions are influenced mostly by frequent environmental variations which include 
monthly wind, rain, temperature, food shortage and parisitism. 

We believe that the ideas and methods used in the present article can be 
applied to other nonlinear random systems in biology and physics. 
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